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Ee ENT RODUCTION 


With the advent of the computer, there have been 
tremendous advances in the field of numerical weather 
prediction. Although there are several methods available for 
solving the prediction equations, the main thrust of the 
development of the numerical models has been with the finite 
difference method. Parallel developments in the engineering 
fields have used other methods as well as finite difference 
schemes. One recent technique being used is the finite 
element method. The purpose of this paper was to develop a 
finite element application to a barotropic primitive equation 
model. The objective was to learn the characteristics of a 
finite element model and compare it with a similar finite 
difference model. eee initial conditions were used to 


insure the most accurate analysis possible and to simplify 


the comparisons. 


PS 


Pies EN ITE ELEMENTS 


Although the finite element method has only recently 
emerged as an efficient means of solving differential 
equations, its beginnings can be traced loosely back to early 
times. The basic concept of the finite element method is 
that a solution can be accurately determined by a sum of 
Simple, easily computable functions. This procedure is not 
new. Martin (1973) states that a Chinese engineer in 
A. D. 480 was able to determine upper and lower bounds on 
the value of mt to seven digits. This was accomplished by 
accurately determining the area of a circle with slender 
inscribed and circumscribed polygons. From calculus of 
variations, the classical Rayleigh-Ritz method shows how to 
approximate a solution to certain problems with a linear 
combination of any set of linearly independent functions. 

The Rayleigh-Ritz method determines the weights that are 
associated with each function. For a great many years both 
engineers and mathematicians were hung up on the idea that 
each of these functions must be essentially non-zero over the 
entire domain of the problem. For a large complex solution, 
the determination of the weights could be prohibitive. Then 
in the early 1950's engineers applied the Rayleigh-Ritz method 
Mpedividing the entire domain into many smaller pieces or 
elements. Hence the anachronism, finite elements. In this 


respect, the spectral method and the finite element method 
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are guite similar. The spectral method is a combination 
of sines and cosines defined over the entire domain while 
the finite element method is a combination of low order 
polynomials, each polynomial or basis function, being non- 
zero only over a small "finite element." The question is 
then how to choose the coefficients for the polynomial to 
best approximate the answer to the equation. The classical 
Rayleigh-Ritz procedure for linear self-adjoint problems 
formulates the problem as the solution to a minimization 
of apositive definite functional. The coefficients then are 
chosen to minimize the error in using a finite number of terms. 
Galerkin proved that the same coefficients result if the 
problem is formulated as a self-adjoint linear differential 
equation and then the coefficients are chosen in the following 
manner, Insert the linear combination of N basis functions 
into the equation, multiply this resulting equation by N 
linearly independent test functions (usually the basis 
functions themselves) and integrate to get N equations for 
the N coefficients. (Observe that if the test and basis 
functions are the same, multiplying the basis and test functions 
involves the integral of the squares of the functions, and 
thus this procedure is related to the least square error — a 
fact which Galerkin proved.) 
One may suspect that this same procedure would also give 
a reasonable error for nonlinear, non-self adjoint problems. 
In summary, the Rayleigh-Ritz-Galerkin procedure for 


nonlinear, non-self adjoint equations involves subdividing 


1S 





the domain into a set of elements, approximating the dependent 
variables by a linear combination of low order polynomials 
(having value zero except over the particular element), 
inserting these approximations into the equation, multiplying 
the equation by a test function (also a low order polynomial 
having a one-to-one correspondence with the basis functions) 
whose purpose is to minimize the residual between the 
approximate solution and the actual solution, integrating 
over the entire domain, and finally solving the system of 
equations assembled during the integration for the solution. 
Since the basis functions are zero over almost all of the 
area over which integration takes place, this procedure 
produces a system of equations for which the matrix is very 
Sparse. The Rayleigh-Ritz-Galerkin method has become the 
most popular finite element method and is used in this paper. 
emeean (1973), Zienkiewicz (1971), Strang (1973), Schultz 
[myo ), Norrie (1973), and Aziz (1972) give in-depth 
theoretical descriptions of the finite element method. 

The first step in the Rayleigh-Ritz-Galerkin method 
requires the division of the domain into a set of elements. 
The element shape in this paper was chosen to be triangular. 
Section III contains a description of the subdivisions 
mtcilized. 

The next step requires representing the dependent variables 
as a linear combination of low order polynomials. In this 
paper the polynomials used for both the test and basis 


mci lons were linear in A and 8. A linear basis function 
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over a particular element can be visualized as a plane with 
value one at one of the vertices of the triangle, value zero 
at the other two vertices of the triangle and identically 
zero over the rest of the domain. Since any one element has 
three points, there are three planes used in approximating 

a function over each element. This can be seen in Figure l. 
Any function, say u, can be represented over the particular 
element as a linear combination of the three planes shown 


in Figure 1 by the equation 


3 
Velement — a a (t) ie Geet) 
j=l 
where 0. represents the scalar value of u at point j of the 
element and a are the basis functions. At the boundary of 
two adjacent elements, the dependent variables are continuous. 
It may be noted, however, that since only linear functions 
are used for the representation, derivatives are not 
continuous along the boundaries. 

After the approximation for the variables have been 
substituted into the equations, the next step involves 
multiplying by a test function. Since in this paper, the 
test and basis functions were the same, Figure 1 also shows 
the three planes making up the test functions over any one 
element. 

The next step in the Rayleigh-Ritz-Galerkin procedures 
is to integrate over the domain. Since each basis and test 


function is zero over the domain except over the particular 


Ny 




















Three planes representing a function over an 
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element, the global integration can be performed by 
integrating locally over each element. In this manner the 
equations are written and then solved for at each node. The 
basis and test functions may be visualized from the view- 
point of a node instead of an element as shown in Figure 2 
where the six basis functions having value one at point i,j 
for each of the six elements are shown. the S1x planes 
(basis functions) shown make up a "pyramid function" for 
the grid point i,j. This "pyramid" function is the test 
function that multiplies the variable representation as 
described above. If this procedure is repeated for the N 
grid points a system of equations with N equations and N 


unknowns will be generated. 
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(i- 1,j-1) (i+1,j-1) 


PrGUmie2. Pyramid function at point 1,7 


20 








Mita ROLROPTC PRIMITIVE EQUATION MODEL 


The time integration was performed using a Galerkin finite 
element application to the so called shallow water, "primitive" 
equations. A time extrapolated Crank-Nicholson method was 


used during formulation. The Crank-Nicholson method is 


explained in some detail in Appendix A. 


Pee ERIMITIVE EQUATIONS 
The primitive equations, in spherical coordinates, for 


this model are as follows: 


80 4 — 2 (gu) + say & WV cos 8) = 0 (III=1) 


du u Ue 2 OU uv a : 
Meancos 6 31° a6. a 2 9 - 28 sin @ v 
Sieg, es cK) = 0 (TII-2)} 
a cos 6 3A 
oV u ov VY ov u? 
st + Roos ob ay t ape tm tan @ + 22 sin 6 u 
+ L 2¢ = 0 (ITII-3) 
a dé 


After equations (III-1l), (III-2), and (III-3) were 
expressed in finite element form, they were solved in the 
Bemlowing Order: the continuity equation (III-l), the 


u-egquation (III-2), and the v-equation (III-3). Each equation 


ai 














was solved by an iteration procedure until a converged 


solution was achieved before proceeding to the next equation. 


B. GRIDS 

Two different grids were used for experiments during 
development of the computer model, namely, the simple 4,86 grid 
and the icosahedral grid. The iA,®@ grid conserves angular 
distance in radians while the icosahedral grid nearly 


menserves linear distance in meters. 


fee x,o Grid 

Since the primitive equations are normally expressed 
in the spherical coordinate form, it seemed natural to use 
a grid simply subdivided by longitude and latitude. 
Longitude ranged from zero to 27 and latitude, from - 5 to 5 : 
The grid was then evenly subdivided into the desired intervals. 
A ten-degree interval generated a "Square" with a side of 
length ten degrees. The construction of one diagonal across 
each square further subdivided each square into two triangles 
with the same square radian area. (See Figure 3.) The 
intersection of adjacent triangles' apices defined a node. 
Each node was then assigned a global correspondence number. 
The numbering scheme ran along constant latitude lines in 
order to keep the bandwidth of the coefficient matrix as 
Narrow as possible. Cyclic continuity was maintained by 
numbering the nodes at A = 217, with the same number as the 


4 = 0 node given the same latitude circle. Evéry internal 
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node was supported by six soe nodes; the north and south 
pole nodes were supported by four other nodes. A drawback 
to this grid was the north and south poles being represented 
by a line rather than a node. The primitive equations in 
the spherical coordinate form are singular at the poles. 
Thus this grid placed many nodes at one point (the north or 
south pole) where the equations are singular. 
2. Icosahedral Grid 

Williamson (1968) and Sadourny, et al. (1967) 
pointed out the advantages of a grid which is nearly 
homogeneous with respect to area, and both described methods 
of generating a grid which used an icosahedral spherical 
figure. Cullen (1974) also used a regular icosahedron while 
integrating the primitive equations. 

A regular icosahedron on a sphere consists of twenty 
major spherical triangles with twelve vertices (see Figure 4). 
An icosahedral grid is then superimposed by subdividing the 
major triangles into smaller triangles of nearly uniform 
area, which is the most important feature of this grid. If 
a model can be run on a global grid with nearly equal area 
triangles then the problem of decreasing time steps 
associated with decreasing Ax as the poles are approached 
can be eliminated. Cullen (1974) subdivided each major 
Spherical triangle along latitude/longitude circles but 


pointed out that a great circle subdivision would best 


conserve equal areas. 
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FIGURE 4. Regular icosahedron 











The icosahedral grid used by this model consisted 
of a subdivision of each of the twenty major spherical 
triangles. Each major spherical triangle side was divided 
into n segments and the points were joined by arcs of great 
circles to produce a grid as shown in Figure 5. The 
mathematics are shown in greater detail in Appendix B. The 
number of points, N, in the grid where each side of a major 
spherical triangle is divided into n pieces is given by the 


following formula 


Nae — (2k0n= + 2 (IITI-4) 


The number of triangles, T, generated by the same grid is 


given by 


The Me 20n? (II1I-5) 


With the exception of the vertices of the major spherical 
triangles,each node was supported by six other nodes. At the 
twelve major vertices, the nodes were supported by five nodes. 
The global correspondence number given each node was assigned 
by starting.at the north pole and going around the latitude 
band as shown in Figure 6. The points of the triangles lying 
On zero and 360 degrees longitude were assigned the same 
global correspondence number to maintain cyclic continuity. 
This can be seen in Figure 7. 
3. Global Correspondence Table 
A Galerkin finite element approach requires inner 


products of the dependent variables, say f, with a pyramid 
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FIGURE 5. Icosahedral grid 

















NORTH POLE 





SOUR FOLe 


FIGURE 7. Elements displaying cyclic continuity 














manction, V,- The inner product in spherical coordinates 


%#S defined as 


<£(A,8),V,> = SS £(A,8) V, a* cos 8 ddde (aia) 
global 

fies pyramid function Vie at any point i, has value one at i, 
linearly decreases to zero at the surrounding points and 
remains zero over the remainder of the globe. The global 
integration can be performed by integrating the particular 
function on each triangle, and placing the value of the 
integral in the proper place in the global matrix. A further 
explanation of the interaction between basis functions and 
test functions can be found in Section III.C.1. The global 
correspondence table is a means of properly scattering each 
triangle's surface integral into the appropriate space in 
the global coefficient matrix. 

Given a triangle with local coordinates 1, 2, and 3, 
the inner product over the triangle will produce a value for 
the three points 1, 2, and 3. Point 1 will receive the 
values from the interplay with itself (1,1), with point 2 
(1,2) and with point 3 (1,3). Point 2 will receive values 
meemernterplay with (2,1), (2,2) and (2,3). Similarly point 
[ewe receive values from (3,1), (3,2) and (3,3). Each 
triangle has a local 3 x3 matrix which needs to be scattered 
meer its global position in the coefficient matrix. Since 


each of the N points are multiplied by N global pyramid 
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functions, an N x N global coefficient matrix is formed 
during a complete integration over the globe for all points. 
Each row, 1, in the matrix represents the equation for point 
i. A table had to be set up which correlates the local 3 x 3 
Matrix to its place in the global N x N matrix. The first 
step in developing a global correspondence table is to number 
all the nodes. The numbering scheme should be chosen so as 
to minimize the bandwidth of the coefficient matrix. For 
this paper, another approach was taken which reduced the 
storage requirement to a Minimum (Section III.C.3). George 
(1971) includes a subroutine which, using any grid, will 
optimumly arrange the correspondence table to provide minimum 
bandwidth. 

The icosahedral grid is indexed by starting at the 
north pole, point 1, and numbering around the next consecutive 
latitude band until all latitude bands are completed. The 
next step is to impose a local 1-2-3 coordinate onto each 
triangle and compile a local versus global table. Table I 
shows the first 10 triangles and its global correspondence 
table. The correspondence table is kept in a matrix and 
called during area integrations to scatter the 3 x 3 local 
matrix into the proper place in the coefficient matrix. The 
Computer Program section contains a program to generate the 


correspondence table. 


on 








TRIANGLE LOCAL COORDINATE 
2 3 


NUMBER 





TABLE I. Global correspondence numbers for 
the first ten triangles 
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Seer iLNITE ELEMENT APPLICATION 
ma oreasintegrations 

The inner product Wa Nae can be computed by performing 
the necessary integration over each triangle separately and 
Meeecrlouting the values to the proper position in the matrix. 
The actual integration can be done by a numerical scheme such as 
Gauss quadrature, or by use of an analytic approach. Due to 
the low order polynomials (linear) used, exact expressions 
may be derived for the integration of any function over an 
element. These formulas are used through a coordinate system 
called area coordinates. Strang and Fix (1973) state that 
area coordinates are known to engineers as triangular 
coordinates and to mathematicians as barycentric coordinates. 
Desai and Abel (1972) call them natural coordinates. 

In the formulation, it becomes necessary to perform 
the inner product Vag Ua oe Given the three pairs of 
coordinates (x,y), of the triangle in the x,y plane, there are 
three corresponding coordinates (L,,L,,L,) in the natural 
coordinate system. This is shown in Figure 8. Zienkiewicz 
(1971) shows the relationship between a point x,y and the 


area coordinates (L,,L5,L3), as 


x = Lx, + L5xX5 + L3x, CLri= 7.) 
ieee deo * bY 5 ee! 
1 Shae (III-9) 


where Lb, = Meets aD JA; Li = A./A 














and A = total area of triangle; and A, 1A5,A2 = area of the 


smaller triangles. He also shows that 


! 
ff us” ay ey yee OR (XL II-10) 
A (atbt+cat+2)! 


This formula was used to perform the integration. For 


example, the inner product Sepa at j=i=l is 


_ 2 y 2) Or 

UO i JJ Va dxdy = (2+0+0+2)! oh 
eae = is 
= yt 9 GS = 


and the inner product <V50V5> at point j=2, i=l is 


_ wa lb One 

<VoeVy> = JJ VQV, dxdy = —rsreoeayy 2A 
a _ A 
= yy a a a 


With this formula, any necessary product of basis 
functions can be integrated over one triangle. 


If equations (III-7), (III-8), and (III-~9) are now 


solved for Lj Lo, L. and written in matrix form the result 
ZS 
Ly 2A Dy ay 1 
gee _ 
Lo = Si 2A b. a, x (III-11) 
L. 2A De eke Y 
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FIGURE 8. Cartesian coordinates vs. natural coordinates 
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FIGURE 9. Triangle definitions for area coordinates 











as shown in DesSai and Abel (1972) where the a's and b's are as 


defined in Figure 9. Differentiation of (III-11) shows that 


3 ea, 5 
o— = (IITI-12) 
ox aa 2A ob, 
3 te 
Cla ) tt (III-13) 
3Y jan 
These two equations will be used when evaluating the 
derivatives of the basis and test functions. 
For example, assuming V5 = ib , the derivative Vay 


meepoint J] is 


et Ge ae Pas ree | 
1x 2A OL, 2A OL, OA aL. 
but 
dV dV 
my 7° + pet 
2 3 
OV, 


and i = 1 at point 1. Conseguently, 
1 


Oo 


_L 
1% IR 


The inner product SU at point j=2, i=l becomes 


bo 
=~ V, dxdy 


< = 
cee 2A ‘2 


2x" 1 


Piso 0! 
(1+0+0+2)! 2h 


Aes 


OF 
fh 
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eas with these simple formulas, all the inner products can 
readily be evaluated. 
2. Nonlinear Terms 
The nonlinear terms in the system of equations were 
quasi-linearized by the time extrapolated Crank-Nicholson 
method. When a uu, term was encountered, it was replaced 


by u*u, where 


(ee Jee 


Nf 
G 
l 
Nl eH 
G 


Thus when solving for time level N+l, all the * quantities 
will be known. In this manner the nonlinear terms are quasi- 
lienarized in time to facilitate integration. The inner 


r 


this paper. In matrix notation this becomes 


* 
peodauct ee a represents u*u, in the notation of 


= * - 
OA, OSE V Va VG? (TiE= E25) 
mmethe local 3 x 3 matrix, any point within the local matrix 
becomes 


= <q*V_V. Sie (IITI-16) 


; .> + <qk : .> + <qt* , 
4 11 peal C2 VaV V5 Oe Va 


DO or BN Ge 
In this manner the known * functions are integrated into the 
Galerkin space. The nonlinear terms are "linearized." Each 
equation becomes one equation with one unknown. The three 


equations are coupled when the new information at time level 
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N+l is used to solve the equations at time level N+2. A 
constraint on a linearization scheme is that a large time 
step cannot be taken. The time step should not be so large 
aS to cause the change in the variable to be larger than the 
feemcation error of the approximation. This linearization, 
which uncouples the three equations, causes some inconsistency 
with respect to the three dependent variables. 

3. Matrix Storage 

During a global integration, an N x N coefficient 
Matrix is generated which is very sparse due to the fact that 
the maximum number of triangles supporting any one point is 
Six. These six triangles provide interaction between the 
seven points involved. Each row, i, in the N x N matrix 
therefore represents the equations written down at point i, 
and each row would have, at a maximum, seven entries. If 
the matrix were condensed by omitting all terms identically 
zero, the size would be N x 7. This represents a sizable 
reduction in core storage. 

The method of condensed matrix storage offers the 
advantage of minimum core storage, For large fields, a 
complete N x N matrix can easily exceed core capabilities 
of most computers. The storage of a sparse matrix would be 
extremely wasteful, hence, the condensed method was adopted. 
If the coefficient matrix were symmetric or had a narrow 
bandwidth, then storage could be easily accomplished using 
a smaller amount of computer core. The icosahedral grid 


offers neither a symmetric matrix nor one with a narrow 
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bandwidth. The disadvantage of the matrix storage scheme 
was the need to search a correlation matrix when a value 
needed to be placed in the coefficient matrix. An efficient 
searching routine was developed and can be seen in the 
Computer Program section under Subroutine SEARCH. But, 
unfortunately, this subroutine had to be executed many times 
during the time integrations. 

PimeEmeoemodel, an iterative procedure was used to 
solve the coefficient matrix. In order to reduce the N x N 
mecrix toan N x 7 matrix, a correlation table had to be 
developed. The correlation table, a separate N x 7 matrix, 
was assembled which contained the seven points involved in 
any one row of the coefficient matrix. An example will show 
more clearly the process involved. Table II and Figure 10 
show the triangle, points and correlation table involved with 
point 2. In an N x N matrix, point 2's equation (row 2) would 
meme nNOn-zZero entries in position (2,1), (2,2), (2,3), (2,6), 
fy), (2,8), and (2,16). For example,in the N x 7 coeffi- 
Senne Matrix, entry (2,16) would be stored in (2,7) of the 
condensed matrix and entry (2,8) would be stored in (2,6). 
Bernd Matrix multiplication, the proper address of a vector 
component had to be selected as well as the proper address 
in the N x 7 matrix. This can be accomplished by calling, 
mee example, position (2,7) of the correlation matrix to 
find that entry (2,7) of the coefficient matrix goes with 


entry 16 of the vector. The process is reversed when it 








Legend Q) Triangle number 1 


1 Point number 1 





Figure 10. Diagram of triangles l, 5, 6, 7, 19, 20 and 
modes el, 2, 3, 6), 1,785 16, 


PABEE Tl. Correlation table for nodes 2, 7, 8 


CORRELATION TABLE 


N x 7 Matrix 
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Mmecomes necessary to scatter the local 3 x 3 matrix built 
during each triangle's surface integration. The global 
number of each point in the triangle determines which row 
into the N x 7 coefficient matrix the value will go. The 
column is found by doing a search of the row from the 
correlation matrix until the second number is found. The 
column of the correlation table in which the second number 
was found determines the column in the N x 7 coefficient 
Matrix. For example, triangle 6 will have a local 3 x 3 


Matrix with global numbers which looks like 


(2,2) (2,7) (2,8) 
(7,2) (7,7) (7,8) 
(8,2) (8,7) (8,8) 


Row 2 of the coefficient matrix will contain the values of 
meres, (2,7), and (2,8). A search of correlation matrix's 
row 2 shows that (2,8) of the N x N matrix goes into (2,6) of 
the N x 7 matrix. 

The correlation matrix is developed by searching the 
global correspondence table for the six triangles that contain 
point 1. These six triangles are then compared and sorted to 
produce the seven points interacting with point i. The 
seven points are then arranged in ascending order. 

A copy of the program which produced the correlation 


table can be found in the Computer Program section. 
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IV. INITIAL CONDITIONS 


The initial conditions used in this paper were computed 
from an analytic solution to the nonlinear balance equation. 

Haurwitz (1940) developed a stream function and showed 
that harmonic waves computed from this stream function will 
move with constant angular velocity in a non-divergent baro- 
EEO DLC atmosphere. The stream Eaneeent w, used by Haurwitz 


is given by 
y = A* sin(m) - vt) sin 6 cos m@ - Ba’sin 6 (IV-1) 


where A* and B are constants, mis the wave number, v is the 
angular wave velocity, a is the radius of the earth. 
The constant B is related to the wave number and angular 


wave velocity by 


ee acs) ee ee Ty - 
= Seba CELE a: 


si< 


where (v/m) is the angular phase speed, 2 is the angular 
velocity of the earth and M = m+tl. To allow a phase speed 
propagation comparison with a finite difference scheme, the 
constant A* was arbitrarily chosen to be the same as the A 
mea M.S. thesis by Monaco (1975). Phillips (1959) used the 
Haurwitz stream function, wy, in the nonlinear balance equation 


to determine 6. The nonlinear balance equation is 
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2 2 Ve 
Wo = £V-y - Vurve - 2(2-8 2 + 2( 9 ¥ OY Hees 
dX0Y 2 2 
dx OY 
which in spherical coordinates becomes 
il iL an i) ,) 3 ¢ 
Mm 4.2, 4,2 Be SC eoee 
a cos” 6 dA cos 8 398 
f Le BAG ilpenn av 
° >) ae ae pereiede y a 
a cos” 6 dX cos 8 38 
1 oy a¢£ len Couey 
> Se 2 ( ) ———-) 
a 086 080 a’ cos 8 3gAd8 
2 2 
2 0 pow 2 
+ 5 =f (IV-4) 


Phillips (1959) found the solution to the nonlinear balance 


equation for the geopotential distribution, ¢$, to be 


o = a? Ee) ct a? B(6) sin(mdA) + a“c(6) (2 sin’ mA - 1) 
(IV-5) 
where 
2m 
* 
A(6@) = 2(20+B) cos“6 + 7(=5) “cos 6 aml) eos’ 6 + (ome onean 
a 2 
Z —~n |] (iy Seal 
cos 86 
* 
Digeis) sae 
B(6@) = a’ So [ (m?+2m+2) ~ Gavel ~ aoe 6 j (IV-5.2) 
~  (m+1) (m+2) ae : 
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e(40)? Soa Q [(m+1) cos*9 - (m+2) ] CEN 03) 


a 


and ee) 


The foregoing geopotential disturbance, », must be added to 
the mean height for the level desired. The mean height was 
obtained from the NACA standard atmosphere (Haltiner and Martin, 
1957). Although an attempt has been made to apply the model 
at the level of non-divergence, the equations used in the model 
permit divergence. Phillips (1959) found that the presence of 
divergence will cause the waves in the height fields to move 
slightly slower than a non-divergent model, especially for 
small wave numbers. 

Analytic initial winds can be obtained from the stream 


function, w, as follows 


u= - a (IV-6) 


Ue 


and 
_ JS ee Z 
Vv = = Sees Oi (IV-7) 


From wW, in IV-l1, the values of u and v are 


u=- = [A*s in (mA- vt) costttg ~ mA*sin (m\-vt) cogir lt, stace ~ Ba‘ cos 6] 
(IV-8) 
v= = [A*m sin 9 cos Le cos (mA-vt) J (IV-9) 


These analytically determined winds from the nonlinear 
balance equation were chosen because they are in approximate 
gradient balance and hence more realistic than simply 
geostrophic winds. Subroutine SOLUT of the main program, 
Gemputer Program section, produced the initial fields for 


the latitude and longitude of each point. 
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V. WAVE ANALYSIS METHOD 


A harmonic analysis was performed around constant 
latitude circles to determine the phase speed and amplitude 
of meteorological waves. The initial fields were analyzed 
first to provide a reference condition. 


The Fourier series used was 


F(x) = A. + ) (A cos mx + Bein mx) (V-1) 


which can be represented as one trigonometric function given 


by 
a) = Co + 2 C_cos (mx-6_) (V-2) 
where 
ea vi 
Cn ~ Sin 6 ~ Cos 6 (V~ 3) 
m m 
and 
-J] an 
= = V- 4 
a tan A ( ) 


In order to analyze phase speed propagation, nonlinear 
instability and energy conservation, numerous wave numbers, 


amplitudes and phase speeds were examined. 
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Wiese CRIED VCONVERS TONS 


The grid that is most compatible to this finite element 
model is the icosahedral grid. But this grid is incompatible 
with the grids needed to perform harmonic analysis or 
graphical display. This fact could be a serious drawback to 
an operational finite element model. An advantageous 
characteristic of the finite element eho’ is that the 
solutions are continuous over the element. For this model, 
linear functions were used. To convert the icosahedral 
solutions to solutions at different points involves evaluating 
the function over the element at the point desired. This is 
where the finite element method will have an advantage over 
the finite difference schemes. Finite difference schemes 
must interpolate from discrete point values to an interior 
point. In this finite element method, the basis functions 


are of the form 
f = a+ bA + e8 (VI-1) 


The value of the funciton is known at three points. These 


three equations can be written for points 1, 2, and 3 


fy = sha: bi, + e8, (VI-2) 

f. =d+t bi. + eb. (VI-3) 

f. =d+t bh. 4. eo. (VI-4) 
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We have a system of three equations and three unknowns. The 
value of the coefficients d, b, and e over any one element 
can be easily found. Once the coefficients are known, the 
value of the function at any point A,0 in the element can be 
found by evaluating the polynomial given in equation (VI-1l). 
For this model a 72 x 35 regular 5° longitude-latitude 
grid was required for harmonic analysis and graphical output. 
A subroutine was developed which correlated each point in the 
72 x 35 grid to its appropriate triangle in the icosahedral 
grid. This is subroutine SURVEY in the main program in the 
Computer Program section. The value of the coefficients d, 
b, and e in equation (VI-1l1) were found for each triangle by 
subroutine EVAL in the main program in the Computer Program 
Section. The value of the function at any point in the 
72 x 35 grid was evaluated by subroutine DISPL in the main 


program in the Computer Program section. 
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Vir EXPERIMEN Es 


Experiment 1. The model was run with the same conditions 
using the icosahedral grid and the \,@ grid. The purpose 
of this experiment was to compare a varying Ax length grid, 
Memmost nOn projection finite difference models use, with a 
near constant Ax grid such as the icosahedral grid. The 
conditions used were a 10-minute time step, wave number 4, 
amplitude times wave number = 28. x 10’ meters and a phase 
speed of 10° longitude/day. 

Experiment 2. The model was run on the icosahedral grid 
with wave numbers 4, 8, and 12. The phase speeds were all 
10° longitude/day. A 10-minute time step was used for all 
runs. The amplitude of the waves was changed so that the 
product of amplitude and wave number remained 28 x no) meters. 
This will constrain the order of magnitude of the north-south 
component to remain the same. The purpose of this experiment 
was to observe phase speed propagation of the various waves 
and compare them to the phase speed propagation characteris- 


tics of various finite difference schemes. 


Experiment 3. The model was run with increasing time 
steps with wave number 4. The other conditions were the same 


as in experiment number 1. Time steps of 10, 12, 15, and 


18 minutes were run for a full 72 hour prognosis. 
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Vinee 


PEePpeceimemerl. Initially the model ran on the A,§ grid 
but after twelve hours the harmonic analysis near the polar 
regions showed increasing amplitudes in the high wave numbers. 
Eventually the solutions in the polar regions caused the 
model to become unstable. There are several possible 
explanations for this instability. From a stability criterion 
viewpoint, the Ax near the poles was so small using the i, 84 
grid that the 10-minute time step exceeded the stability cut 
Sere point. This is another major drawback of this grid. 
The decreasing Ax as the pole iS approached is a common source 
of problems for normal finite difference schemes. From 
theoretical considerations, the finite element method should 
be well suited to an icosahedral grid and this should 
alleviate the converging Ax problem. The iX,@ grid has 36 
points around each latitude circle. It can therefore resolve 
wave number 18. Any truncation error created while solving 
the matrix equations will appear as high frequency noise. 

Another reason is that the equations are singular at the 
poles. In order to avoid these singularities, the model does 
fet predict u, v or @ at the poles. Instead u and v are set 
to zero initially and are kept zero throughout the entire 
integration. A value for ¢ at the poles is found by averaging 
the first latitude circle down from the pole and then setting 


the pole and the next latitude circle equal to the average 











value. This will flatten the cap over the pole and make 

the u and v assumptions at the poles more consistent, although 
Obviously this treatment of the poles is not completely 
desirable. 

Experiment 2. Figures 11, 12, and 13 show the results 
of the phase speed propagation test for wave numbers 4, 8, 12 
respectively. While identical tests using finite difference 
schemes were not available, Maher (1974), ina Master's TheSis, 
performed phase speed propagation analysis with similar wave 
numbers and amplitudes using second and fourth order finite 
difference schemes with a comparable number of points. The 
comparison showed that the finite element method was more 
@ecurate than second or fourth order differencing for the 
cases examined. 

Due to the coefficient chosen, the amplitude was minimal 
in the waves at higher aeenee On the right side of each 
figure iS a value showing the number of points around each 
latitude band. As the number of points per each latitude 
band decreases, the phase speed propagation decreases. The 
finite element method therefore has the same relationship as 
finite differencing where phase propagation 1s concerned. It 
takes more points per wave to improve the phase propagation. 
Charts A-I show the wave propagation in 12 hour increments 
for the three wave numbers tested. 

Experiment 3. With the three equations being treated 
separately at each time level and the main forcing functions, 


coriolis force, and pressure gradient force, being treated 


Sl 





explicitly, it is reasonable to expect a linear stability 


criterion approximately of the form 


mi? 


Q 
miles 
“cr 

[A 


where c = the phase speed of the fastest waves and Ax the 
Minimum Ax in the grid. The actual stability criterion for 
this model on the icosahedral grid has not rigorously been 
worked out. However, in notes by Archer (1975), the stability 
Criterion for a simple wave equation was worked out and found 
to be comparable to that with similar finite difference 
schemes. The assumed stability criterion predicted a maximum 
time step of 19 minutes. It was found that the model with 

an 18 minute time step became unstable after 12 hours and 
after 24 hours with a 15 minute time step. However the model 
remained stable with a 12 minute time step up to 48 hours. 

In all cases, instability in the equatorial regions was 
observed after a sufficient period of integration. For the 
wave number 4 case, this instability occurred after 84 hours 
of real time integration. The instability occurred sooner 


for larger time steps. 
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FIGURE ll. 


Phase angle (degrees longitude) vs. latitude 
for icosahedral grid, wave numb¢x 4, phase 
speed 10°/day and A* = 7.0 x 10°. (Latitudes 
with near zero wave amplitude are not included 
and time is given in hours.) 
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Chart A. Initial PHI field analysis, wave number 4, phase 
speed 10°/day, A* = 7.0 x 1022 
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Chart B. 12-hour PHI field forecast, wave number 4, phase 
speed 10°/day, A* = 7.0 x 107. 
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Phase angle (degrees longitude) vs. latitude for 
icosahedral grid, wave number 8, phase speed 
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BEGURE 12. 


10°/day and A* = 3.5 x 107. 


(Latitudes with 


near zero wave amplitude are not included and 


fame LS given in hours.) 
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Saart D. Initial PHI field analysis, wave number 8, phase 
Speed 10 /day, A* = 3.5 x 10’. 
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Mart &. L2-hour PHI field forecast, wave number 8, phase 
speed 10°/day, A* = 3.5 x 10°. 
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Chart F. 24-hour PHI field forecast, wave number 8, phase 
speed 10°/day, A* = 3.5 x none 
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FIGURE 13. 


Phase angle (degrees longitude) vs. latitude 
for icosahedral grid, wave number 12, phase 
Epc /dayeand f° 9" 2a. l0’. (Latitudes 
with near zero wave amplitude are not included 
and time iS given in hours.) 
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Chart G. Initial PHI field analysis, yave number 12, phase 
speed 10°/day, A* = 2.3 x 10°. 
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Chart I. 24-hour PHI field forecast, wave number 12, phase 
speed 10°/day, A* = 2.3 x 107. 
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IX. CONCLUSIONS 


This model showed that the finite element method is 
applicable to a meteorological system and is competitive with 
Finite difference schemes. The size of the time step that 
can be taken by a finite element method is comparable with 
the one for a finite difference scheme. The accuracy of the 
phase propagation of the finite element method was closer to 
Bic analytic solution than second and fourth order finite 
differencing for the comparisons made. 

Based on the increased accuracy and comparable time steps, 
future research should be continued on the finite element 
method. Specifically, more advanced methods of solving the 
system of equations are available and should be utilized, for 
example, ADI or a cyclic reduction method. Also higher order 
polynomials should be used for the basis functions. In fact, 
the real advantage of finite elements is not realized until 
higher order polynomials are used, from which greater accuracy 
can be expected. 

Research should continue in the polar regions to ascertain 
the ability of the method to resolve flow over the poles. 

This will clearly improve the model as it will have realistic 
physical conditions at the poles instead of artificial 
boundary conditions. To extend the allowable time step, 

the equations should be written so that the three equations 


would be coupled at any one time step. Also the instability 





observed in the equatorial regions after long time integra- 
tions should be further investigated. 

The barotropic model should be expanded into a multi- 
level baroclinic model with a heating package. Real data 
should also be applied to the barotropic and baroclinic 
models to test the method's ability to handle actual 
Genaitions. 

Finally, research should be performed using smaller 
elements perhaps selectively in certain areas. With the 
Success with the icosahedral grid, the application of the 
finite element method to fine meshed models shows promise 


to provide greater resolution in specific regions. 
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APPENDIX A 


EQUATION FORMULATION 


In the initial stages of this paper several decisions 
were made as to the order and method of solving the three 
coupled equations. It was decided to uncouple each equation 
so that it would be one equation with one unknown. Each 
equation would be written in the time extrapolated Crank- 
Nicholson method. Although this method has a slower 
theoretical convergence rate than solving the three equations 
coupled together, the number of computations per time step is 
reduced. If the three equations are coupled, then each 
equation's solution after one iteration at one time step is 
used in solving the other two equations for the same iteration 
for the same time step. This requires area integration of 
each equation for each iteration at any one time step. For 
example, with the Crank-Nicholson method the ¢ in the u- 
equation is written at time Nts . Since this is an unknown, 
each iteration of the $¢ equation for time, Ntl, gives a 
different te: . The pressure gradient term in the u-equation 
would have to be integrated again as it is not the same as 
it was at the last iteration. The uncoupling of the three 
equations was chosen because it sharply reduces the number 
Of computations per time step. The advantage of coupling the 
three equations at any one time step would be that the 


equations would be more accurate and consistent and larger 
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time steps would be possible. Even with the restriction on 
the size of the time step, a 12-minute time step was possible 
before instability was experienced. 

The first equation solved during one time step is the 


continuity equation 


16 ee se eo Z - 
ane * 2 cos 6 OA (gu) + a cos 6 96 (gv cos 8) 0 (A-1) 
An inner product is now performed on each term with a 


global pyramid function Vis The inner product is defined 


<E(X,8),V,> = ff £(d,6) V, a” cos 6 dddg (A-2) 
global 
The at is not carried since it 1S a common factor and can be 
cancelled now. 
The continuity equation becomes 
ia aL AL 3 


Cos 6b On (ou) pa a5 a aC ceral 5 76 (ov cos 6) Ee 





= QO C3) 
If the second and third terms are integrated by parts, 
the continuity equation becomes 


Zan ov. 
Ee ee ee LT 





2 
ad a’ 
< e V5 > + = if (ouv, ) 








0 
T 
OV. 
a 2 1 $v cos 6 aL 
fy = = 
= lit Om COS 8 7) ar ais =< REE > 0 
g 2 (A-4) 
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The second term in the above equation is zero due to 


meric Continuity and the fourth term is zero due to the value 


O£ cos 6 at - > and > - The continuity equation becomes 
OV. OV. 
ad aoe ou ee LL ea ores, <3) a a = 
Pat’ “i? a “cos 0’ DN a “cos 0’ 30 ~ : ee) 


To extrapolate the dependent variable forward in time the 
Crank~-Nicholson method first uses a simple forward difference 
in time to N+l followed by an average at time level N+l with 
values at time level N for the space derivatives. The space 
derivatives, therefore, are evaluated at time level Nt. 
Since the three equations were uncoupled, any variable that 
is not the variable for the particular equation is evaluated 
at time N+ while the variable for the equation is evaluated 


at times N and N+l. The continuity equation becomes 


N+1 .N N+1 N+ 9V. N N+} 93V. 
oi So ein pees ee u ae 
a At ) ee 2a : GOS om oo) 2a ~ cos ey ON a 
N+1l N+ ov. N N+3; OV. 
iL ae i i: Wen ee 7. 
= a a at ee —a> 
waeecos 0 °°° °' 36° Za “cos 6 ©°S "56 


(A-6) 


The variables u and v at time N+ are represented by the 


second order approximation 


wo 2 => ue ~ = u - > u (A-7) 
vite = ye % = vs 5 ve (A-8) 
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This completes the uncoupling and the continuity equation 





becomes 
N+1 OV. N+1 ov. 
N+l 4N pe d ui 1 d ieee 1. 

< - Vee = = —_ sae a 

(9 me) an 2a ; cos 6 ’d0iX ge Ga Gren 8655 

N OV. N OV. 
_ Ne a1 o v* ah _ 
Deal Es Sa Y “See a Se iy. ° Bop (A-9) 


The continuity equation is now one equation in one unknown. 


Using Galerkin formulation with Einsteinian notation, | 


the variables u, v and ¢$ are represented by 


bl = a. () v, (A-10) 
Vv = 2a ey V, (A-11) 
a ee, ve (A-12) 


where NE are the low order polynomials that are the basis 
functions. In this paper we used linear polynomials of the 


form 


f=d+ ba +eé (A-13) 


The coefficients a, 8 and y are the scalar values of the 


variables and are functions of time only. 


Ty repeated index implies summation with respect to that 
index. 
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Using this new notation, the continuity equation becomes 





* 
N+1 _N Malet Oe ee 
< mee )V.;V.> - — ee See 
CY, Ys) i 2a “Ys cos 0’OA 
* 
a oye cos 8 yy ae eas Siete 4 V4 V4 ra 
J k cos 6 j k’9d6 2a V5 Gos Ur 3 
oV 
a 
n <y38 cos 6 | yy _ i, 295, (A-14) 


k cos 6 j k’ 38 


where Vv; is the pyramid function at the point where the equation 
is being written. It is now clear why the pyramid function was 
introduced. This procedure is repeated for each grid point 
leading to a system of equations for the value of dependent 


at the N grid points. 


variable 14" 
The trigonometric functions in the equation can be handled 
in several ways. Cullen (1974) treated the trigonometric 
functions as constants over each triangle. As the grid length 
decreases the accuracy of this approximation increases. The 
decision was made to interpolate the trigonometric functions 


into the Galerkin space. The trigonometric functions were 


represented by 


(A-15) 


II 
wy 
< 


cos 6 


Salil (3 ag eee (A-16) 


doe 





Ber Case of notation the trigonometric functions were 


cancelled where possible from the inner products and we now 


define a new inner product 


PLY ie 


<u,V;> = u Vi dide (A-17) 


on 


m7 2 


The continuity equation becomes 


7 av 
NHN _ At} Ne1_* 7 
ee) 5% VG? ~ DalSt5  RNVGYR ON 


av aV 
N+1,* i At} oN * i 
es 5a’ 30 =| ~ FalS%5 ¥en TN 
Na * ee 
ers Yigg -] OC ie 


Since everything is known in this equation, with the exception 


ong . and the Y. are functions of time only, we can take 


the Ys outside the inner products and write the following matrix 


equation 
Ntl ig _ Ata) = Nowy Atg : 
Ys [A aa B| = Ys [A + aS B] (A-19) 
where 
a ee yk? Vie. 
and 
a * av. x av. 


dz 





The matrix equation is manipulated to solve for the 


N ; 
change, en eee y 


j 
yl K St 5] = yS(R + == By + VS GE B - ot 5) (A-20.1) 
ie -5£ 5] = Ys(K - st Bp] + Y —= 5 (A-20.2) 
ev IK - c= B] = Y “t B (A-20) 
where 
_ n a A e; (oO) 


The zonal equation was developed in the same manner as 


the continuity equation and the matrix equation derived was 


Ale ce) = Ste = oN Se Ue 
J 2a a J a 
where 
Z pee . A OV. 
= > < > - 
D <a, V, 352 pia SBE, Vy Vy, ee (cele) 
and 
rE <2anN i > < ‘ OV" > 
Pi 2a 0B eV, ge gS eo 4 
<a. Be > 21.2) 
+ OB aN VV VV 52 (A : 
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The meridional equation becomes 


Nex, At eee Nit ieee 
e5 [A es f= 8. ac (A-22) 


where 


< 


. x 9 ; * 


av. 
oul = 
EnYEV, 7g V2 (A= 22745) 


© 


and 
re * & * 
SO Op NYE VE Vie Vi? + K2anoay €5 NV, Vp Vay Va? 
A ea? a aa. 22 
SMe, ayy Sano ae (AS 227) 


The matrices in equations (A-20), (A-21), and (A-22) 
are built by the procedures in Section III.C.1. Once the 
matrices for each equation have been built, the system of 
equations derived can be solved by any conventional process. 
A simple Gauss-Seidel iterative procedure was chosen to solve 
the matrix equations with a relative improvement check used 
as a cutoff criterion. Subroutine SOLVE in the Computer 
Program section solved the equations. 

While the purpose of this paper was to develop a finite 
element barotropic primitive equation model, it was hoped 
that a considerable savings in time could be realized. Further 
time savings can be achieved if more sophisticated techniques 


such as, Successive Over Relaxation, SOR, or Alternating 
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Direction Implicit, ADI, are used. Leslie and McAvaney (1973) 
show that the computing time can be greatly reduced by using 
some of the newer direct techniques. These latter were not 
employed during development since programs for these methods 


were not available in the library of programs. 
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APPENDIX B 


ICOSAHEDRAL GRID FORMULATION 


A spherical icosahedron is made up of twenty equilateral 
spherical triangles. Each interior angle is 27/5. The 
method of subdividing the icosahedron is accomplished in 
three steps. First, the top five major spherical triangles 
are partitioned. Second, the Southern Hemisphere triangles 
are reflected from the Northern Hemisphere. The last step is 
the subdivision of the equatorial triangles. 

The first Northern Hemisphere major spherical triangle 
is subdivided and the next four adjacent major spherical 
triangles are given the same solution for longitude and 
latitude of the nodes except that the longitudes are displaced 
the proper multiple of 27/5 radians. In all cases, the purpose 
of the subdivisions is twofold. The main purpose is to specify 
the latitude and longitude of every node. The second purpose 
is to assign a global correspondence number to each node. 

The second step in subdividing the sphere is to mirror 
the Northern Hemisphere spherical triangles into the Southern 
Hemisphere. The solutions are the same except that the 
longitudes are displaced 1/5 radians. 

The third step is the solution of the ten major spherical, 
equatorial triangles. The following description shows how 
one Northern Hemisphere triangle, two equatorial triangles 


and one Southern Hemisphere triangle can equally be subdivided 


76 





and solutions for longitude and latitude of each node be 
obtained (Figure 14). Three laws of spherical trigonometry 


are required. They are: 


Law of Sines 


Sina . sinb 
Sam A» sines (B-1) 


Law of Cosines for Sides 


Bos a = cos b cos c + sin b sin c cos A (B-2) 


Law of Cosines for Angles 


Sos A = -cos B cos C + sin B sin C cos a (B-3) 


where a represents a side (spherical arc) of a triangle and 
A represents the corresponding opposite angle. 

Since all three angles of a triangle are known, the length 
of arc a can be determined by the law of cosines for sides. 
Since ais along a meridian, the length of this arc is the 
Bemacttude of point 2. The longitude of points 1 and 2 is 
arbitrarily chosen as zero radians. Arc ais then divided 
into n segments. The latitudes and longitudes of these points 
can easily be computed. Arc b can be segmented in the same 
manner. Arc ad can be found by using the law of sines with 


angle C and arcs 3-1 and 4-1. Are dis then segmented into 
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NORTH POLE 








the proper number of pieces. The three sides of triangle 
1-4-3 are known, therefore angle D can be found by the law 
of cosines for sides. Arce 1-3, arc e and angle D with the 
iaveor Sines determine the colatitude of point 5, are f. 
Arc 1-3, arc e and arc f£ with the law of cosines for sides 
Perenmane angle E, point 5's longitude. All points in the 
Northern Hemisphere triangle can be solved for longitude 
and latitude by manipulating these laws, arcs and angles. 
All the points in a Southern Hemisphere triangle are 
mirror images of’ those in a Northern Hemisphere triangle 
except that the longitudes are displaced 7/5 radians and the 
latitudes have a minus sign. 
With angle (BtF), arc 1-6, angle C/2, and the law of 
Sines, arc g can be found. This are is then divided into 
n segments. With angle (BtF), are a, are 2-7, and the law 
of cosines for sides, arc 1-7 can be found. This are is the 
colatitude of point 7 but care must be taken as this arc 
can extend across the equator. Arc 1-7, angle (BtF), 
arc (2-7) and the law of sines determine angle (2-1-7), 
mpee Longitude of point 7. Point 8*s longitude and latitude 
can be found by a similar computation using the arcs and 
angles on the other side of the triangle. Arce 7-8 can now 
Mewftound and subdivided. Angle (8-7-1) can aiso be found. 
Then triangle (9-7-1) can be solved to determine the latitude 
made longitude of point 9. The rest of the points in both 
major spherical equatorial triangles can be solved for latitude 


and longitude in a similar manner. The Computer Program Section 


ie, 





contains a program written in Fortran IV for an IBM 360 
computer which will subdivide an icosahedron simply by giving 


it the desired number of segments, n. 
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COMPUTER PROGRAMS 
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